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ABSTRACT 


The feasibility of using a three axis control, momentum bias system with magnetic 
actuators for sun-oriented operations is explored. Relevant equations of motion are 
developed for a sun-oriented coordinate system and control laws are developed for: 
initial spacecraft capture after launch vehicle separation; reorientation from 
Earth-oriented to a sun-oriented operations mode; sun-oriented attitude control; and 
momentum wheel control. 

Simulations demonstrating the stability and time responsiveness of the system are 
performed. Sensor noise input tests are performed to investigate the systems 
susceptibility to imperfect conditions. Cross product of inertia effects are also input to 


test for system instability. 
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DISCLAIMER 


"This thesis contains computer simulation code which has not been tested under all 


possible conditions. Its validity for every possible scenario is not guaranteed." 





Vill 





I. 


II. 


III. 


TABLE OF CONTENTS 


INTRODUCTION... cee debate aes agueiscd baud pooped eaten ya caen eee nesta i] 
Ae BACK GROUND iss srccwaitusecacassceouyeecdia teste tes cones hvaaenpnsaueatinsastenaaaceseeeeaneeesett ] 
Bis SEIS COR Y oi cesepeeticec techies vi ps eee eases eae A ea: py 
PUR POSE sieiicie becca ceaiss cess ten aha esas canes es casts ase espeere eee eee Z 
SPACECRAFT CONFIGURATION AND MISSION ASSUMPTIONS........ 5 
Ac ORB Gcctcictd ievciseteheiaca tector eeeeeeme elena eae weer eI ees 7 
B. BETA ANGLE DEVELOPMENT.............sccsscsssssssssssesceccsssreceessneesessnseesoncoees 7 
C, REFERENCE FRAMEG......:.cccssscccssssssosessscssssscsssnsocessosesensceetsasseesosvenssenneces 12 
1. Earth-oriented Orbital Reference Frame.............csscccsssessereseeeersereeseseeees 12 

2. Sun-oriented Orbital Reference Frame..............ccssecccessecessereerstteeeeseseeeees | 12 

3. Body-fixed Reference Frame............scscssssceneeeersessenseeeseessesensessesesnensenees 13 

D. SPACECRAFT BUS DESCRIPTION. ............:cccssssssessereessssceessstasessesessonsenees 13 
E. ATTITUDE DETERMINATION AND CONTROL COMPONENTS.......... 14 
Vis. “SOMSOTS secs -scozsccsccsides iad canieasienns humucdenaksahin cho ceerssev sereateete pemheea aed sicern eee 14 

@. ‘Barth SemSOr sig.icewscdsssensavesdesstucesdeasnecnspetosdecasoavcvunabiecansguebeestunsessvenisese 14 

by, Sule SCHSOLS seossdehacceviedesccsensscanlsaeadewwcdedacd ie ceaaedsacaes cans cease caneca aeceencsertses 17 

C. Magnetometers...........cscccsccscesssnsssscsessscsssessessersesnsesseessenneasneneeeeees 17 

Di PCUIAt OPS siscecsciseseeiisinscidvasusess dence aveeeave estan sundbuaacanteetaneiveisavepaanapsseeesected ea: 19 

a. Momentum Wheel...............::ccccssssccccsssscesssssscscescssenseceesseesconeseeesonsnseees 19 

b. Magnetic Torque RodS............sccsessessessesecereceaceseessssensssestenssensessensessnens 21 

C. TH rusterrS...............scccssscssssessnsscesssscsesccccessccsessssennsessseecessceccencsacenenscsseess 23 

F. SYSTEM PERFORMANCE ............:::ccccssssseessseessenecsssnceenseseeessseacssnsceesesaeeeees 23 
G. SIMULATION ..........ccccscsessccscesssscccssscecssscescsssecesssceecensaesensnecsensseasonsscescssaareeees 24 
SPACE ENVIRONMENT AND ITS EFFECT ON THE SPACECRAFT........... 25 
A. EARTH'S MAGNETIC FIELD DEVELOPMENT.........:sccssrscsessscessnsseseeceranes 25 
B. TORQUES WHICH AFFECT THE SPACECRAPT.......cccsesesssssssssssrensnenenesens 27 


1X 














E SDistar bance TOPQues ccsiitescvscescesssasiensedcnieeetalecnied sant asudicasenaveenerssle esiaeden 27 


i, PROTO GY MAMIE Cis rsctcsse a) sscsias eet soar beat actacasiaeesenensneectaaleetebtaaeiatac: 28 

Ree, SOAR ESS UMN ai aad ak gees cea enstan ote Sesananansticectusstaaseskeraciacnena: 29 

ommes © 5 028141 81 Si ¢210 0c) 0 Mee re ee ee ne re ee 32 

d. Residual Magnetic Torque............ccccceeeee Seduusesiasceccehasincssantiotaianes 34 

ee: Ota OOS Figen cate ser asst donde a cvnsadessstec tae duadeuieiasan Gavecdeieistern eee a) 

IV. EQUATIONS OF MOTION AND ATTITUDE CONTROL LAWS. ......0 eee 37 
Be BOU RATIONS OF MO TION sviclscsiesssaensicutecreateavaechaelccecedavciwet si evisuneduedesudecbuons 37 

Be “CONT ROD: BOW vp ivgutateaieneet cork casstntet avon aeativcaeauiavtachpieaieaeesieelaiectediccests 39 

C. ATTITUDE CON TROL. DAW 9 iccssaceseapscaccthevtsanecevacdsielnSesebanin eitsaiasedeauaes: 40 

1. Satellite Capture After Booster Separation..............c:ccccccssesscsssscessssesceveneees 40 

2... IWWome nbn WING L-COM OL ognc i oiec ns. veeciadacadacsanaegeastivtsacsibeed itaemovessessen: 44 

3. Momentum Wheel Desaturation Comntrol....0........ccccccsccssssessessscecsessecseseusees 47 

4. Control for Transition from Earth-oriented to Sun-oriented Attitude......... 48 

5;: Sun-oniented Attitude Contd scien ccacsesssansacseesstcesiuiaceisaeeaeaenieeees 49 

V. CONCLUSIONS AND RECOMMENDATIONS. ...........cccccesssccessteeesneesesseeesseeeees 63 
Pre COIN CU US IO N Biers suzesetsnts ota ttt cate nas hi ensenes ts 21a eae tees dete Nao eiatanuehane ee eat 63 

Be RECOMMENDATIONS svccccavecoctece asec advance tires teael dae eeaeceenieeaasa ects: 64 
APPENDIX: SIMULATION BLOCKS AND MATLAB CODE..........ccesssecesseeeess 65 
ETS FOP REPEREN CES io 22ivecen seiarnactcnasauneaesauaveristurervtane eae aetatinenaeaeaceuaens 103 
PINTEVAT OTST RB UTI TS Feige es es reset cre xe est costa et eee eeecae ccna ees 105 





ACKNOWLEDGMENTS 

I would like to thank Professor Terry Alfriend and Professor Harold Titus for the 
time and effort they expended to help me finish this document. Their pursuit of academic 
excellence is commendable. I would also like to thank Professor Leonard of the 
Aeronautics Department for his tolerance of me when I popped into his office for answers 
to difficult questions. The staff of the Servo/Controls Lab at NPS also has my utmost 
thanks for putting up with me for the last three months and reserving me a powerful 
computer (if only they could get Costello out of there). Special thanks goes to Chester 
who felt free to provide his literary expertise every time he got near the computer. 
Finally, I would like to thank my loving wife, Mary, who kept me motivated and put up 


with her grumpy husband who came home late every day. 


Xl 





I, INTRODUCTION 
A. BACKGROUND 

Recently, a National Research Council study group recommended the U. S. 
government develop and launch several satellites to monitor the sun and the effect of solar 
radiation on the Earth's climate. The council's board on global change wants a series of 
orbiters to provide long term, uninterrupted and overlapping observations of the sun. 
Earth-based sensors are unable to detect precise fluctuations in solar emissions because 
the atmosphere absorbs much of the energy. [Ref. 1] 

A sun-watching program spread over many decades may be precluded by a limited 
access to space due to rising costs and limited launch capabilities. Therefore, the group 
advocated the possible use of civilian or military environmental satellites to help survey the 
sun. Satellites under the control of the National Oceanic and Atmospheric Administration 
and the U. S. Air Force Defense Meteorological Satellite Program could possibly be 
brought into service for this purpose, but the study states that satellites dedicated solely 
for this program are needed to avoid interruption in the data being taken. [Ref. 1] 

Candidate satellites for this program would have to have long endurance, a 
relatively high altitude to avoid any atmospheric interference, and a precision sun-pointing 
accuracy of at least one half degree for measuring instrument requirements. A momentum 
bias system with magnetic control actuators appears to be a good, relatively low cost 
choice to satisfy the first and last requirements. 

The popularity of momentum bias spacecraft is due primarily to their ability to 
provide reasonably good attitude accuracy (better than 0.5 degree) with relatively 
lightweight and inexpensive attitude control hardware. Gyroscopic stiffness provided by 
the momentum wheel is a key feature of this system, especially if the momentum bias 
vector can be oriented toward or away from the sun for safe mode recovery. The use of 
magnetic attitude control for satellite applications has become quite popular due to its 
light weight, relatively low power consumption, and long life expectancy. Additionally, 


magnetic attitude control systems provide a smoother torque output as opposed to the 











magnetic attitude control systems provide a smoother torque output as opposed to the 
impulsive output of mass expulsion systems and also avoid catastrophic failure due to the 
lack of moving parts. [Ref. 2] Although magnetic control systems do not produce as 
much torque as mass expulsion or reaction wheel systems, they do provide a viable 
option for satellite applications which do not require rapid attitude reorientations. 
B. HISTORY 

Initial use of passive magnetic control was implemented in 1960 with the Transit 
navigational satellite program. The methodology was to align one of the satellite axes 
along the direction of the Earth's magnetic field. A year later, the TRAAC satellite was 
launched with a system that had the capability of providing magnetic stabilization on 
command. [Ref. 3] Other spacecraft which have used magnetic control include the 
Nimbus, ERTS, ITOS, RAE and OAO satellites. On the Nimbus and ERTS spacecraft, 
the magnets were used to offset fixed imbalances and could be adjusted via ground 
commands. The ITOS system used magnetics for precession control, but the system was 
not autonomous and needed to be controlled from the ground. In contrast, a closed loop 
magnetic system has been used to regulate reaction wheel momentum on the OAO 
spacecraft. [Ref. 2] Finally, the HCCM satellite used a system which is similar to the one 
proposed for use here and the new Iridium system plans to use a momentum bias with 
magnetic control. 
C. PURPOSE 

The purpose of this thesis is to develop, and verify by simulation, the performance 
of an autonomous attitude control system which can serve the needs of the sun-watching 
satellite program. This system should be able to acquire the satellite after separation from 
the launch vehicle, reorient the spacecraft to a sun-pointing attitude via ground command, 
and then perform accurate sun-pointing attitude control throughout the mission life. Most 
attitude control applications that contained magnetic torque rods for control were used in 


an Earth-oriented configuration. This system must be able to perform in an 





-Can past magnetic torque control laws be modified to operate in this new sun- 
oriented configuration? 

-How small can Beta angle become before non-linear assumptions break down 

for the system during Earth to sun satellite reorientation? (Beta angle is defined in 
Chapter II, page 7) 

-Can the desired pointing accuracies be achieved with this system? 


-How long will it take to achieve initial capture and reorientation to the sun? 











Il. SPACECRAFT CONFIGURATION AND MISSION 
ASSUMPTIONS 


The spacecraft bus configuration to be used is extracted from the NPS-STEDI 
project [Ref. 4] developed at the Naval Post Graduate School for the Advanced 
Spacecraft Design class (AA 4871). This provided a baseline for equipment selection, 
disturbance torque estimation and moments of inertia. The development of this thesis is 
not intended to be an extension of the NPS-STEDI project, but a general exploration of 
magnetic attitude control for sun-oriented operations. The NPS-STEDI spacecraft uses a 
momentum biased system to provide three-axis stabilization. Attitude information is 
provided by sun sensors and a single Earth sensor with magnetometer data as a coarse 
backup when Earth sensor data is unavailable. The selected Earth sensor is mounted on 
the momentum wheel, which provides a significant cost saving over using separate 
systems. The system to be developed here will use most of these design criteria, but the 
selection of Earth sensor location will be made after system performance at low Beta 
angles is analyzed. Beta angle is discussed later. 

The type of configuration to be considered leaves the designer a great deal of 
options to explore. The first decision should be which direction to orient the 
momentum-bias vector. It could be oriented perpendicular to the solar ecliptic which 
means that the momentum vector does not need to be reoriented by one degree per day to 
follow the sun. In this configuration, wheel speed would be used to point the solar arrays 
and experiment toward the sun. Unfortunately, this configuration increases the number of 
earth sensors required for pitch angle error sensing as well as the number of 
communication antennae in use. Furthermore, there may be thermal cycling problems 
which arise in this configuration. Finally, the momentum vector would not necessarily 
keep the solar arrays in a sun-pointed direction during safe mode/shut down. Another 
approach would be to orient the momentum-bias vector toward or away from the sun. 


This gives the spacecraft gyroscopic stiffness which would keep the solar arrays and the 








experiment aligned with the sun. This is ideal for safe mode/shut down operation as it 
would help to guarantee power to the satellite after a malfunction by keeping the solar 
arrays pointed in the direction of the sun. However, the momentum vector would have to 
be reoriented by 1 degree per day. Also, unless two Earth sensors are used, pitch angle 
sensing data may not be available throughout the entire orbit period for low Beta angle 
orbit configurations. Ifthe Earth sensor could somehow be mounted on the sun-facing 
side of the satellite, then sensing data would be available during the critical period of 
eclipse (this configuration could cause problems with the mounting of the solar arrays). 
When pitch angle sensing data 1s lost, magnetometer data could be used to provide coarse 
pitch angle attitude sensing for control. 

For this application, the fixed momentum bias vector is oriented along the negative 
pitch axis and 1s used to provide gyroscopic stiffness once the spacecraft is separated from 
the launch vehicle. In addition, the wheel's speed will be used to control pitch angle. The 
momentum wheel size is driven by the need to keep the momentum vector pointed within 
one half degree of the sun line throughout an assumed maximum eclipse period of 35 
minutes. This maximum eclipse time is estimated for a low Beta angle, 1:30 am/pm orbit 
and will decrease as Beta angle increases. 

The next thing to consider in this type of approach, is what type of torque source 
to use for momentum wheel desaturation and roll/yaw control. The two types of systems 
to consider are thrusters and magnetic torque rods. Thrusters will not be used because it 
was decided that no orbit control thruster system was required for this mission, and a 
thruster system would increase weight, complexity, space requirements and decrease the 
satellite's potential life-span. Magnetic torque rods will be used for momentum wheel 
desaturation and roll/yaw control (wheel speed will be maintained at 3000 RPM +/- 10%). 
Roll/yaw corrections will only be accomplished during sunlight periods. Error signals 
generated by sun sensors mounted on the negative Y face (sun pointing) will be fed back 
to the attitude control system to provide torque commands to the magnetic torque rods. 


The torque rods will be used to zero these error signals over the course of the sunlight 


the sunlight periods, thus providing near-zero spacecraft pointing error at the beginning 
of the eclipse period. Individual aspects of this design are discussed below and additional 
assumptions are made. 
A. ORBIT 

The orbit selected was a circular, sun-synchronous orbit with an altitude of 550 
kilometers. This leads to an orbit inclination angle of 97.58 degrees and an orbit period 
of 95.65 minutes. The driving force for this orbit in the NPS-STEDI project was the need 
for an ultraviolet sunrise on every orbit so that a complete atmospheric density profile 
could be measured by the all reflection Michelson interferometer [Ref. 4]. To achieve 
this, the initial longitude of the ascending node was selected to be 127.5 degrees which 
corresponds to a 2:30 a.m./p.m. orbit when right ascension of the sun in the ecliptic plane 
is 90 degrees (summer solstice). For this application, continuous sun coverage is desired 
(i.e. 6:00 orbit) but orbits ranging from 6:00 to 1:30 will be explored to test the validity of 
attitude dynamics assumptions at low Beta angles. 
B. BETA ANGLE DEVELOPMENT 

Before discussing orbital reference frames, the topic of Beta angle and its 
relationship to this problem needs to be discussed. Beta angle is defined as the angle 
between the orbit plane and the sun line, which is the line between the center of the Earth 
and the center of the sun. This angle can be seen in Figure 1 and is defined using Figures 


2 and 3 as 
B = arcsin[SuSeCi + CuSiSw — SuCeSiCw] | (1) 


where S is the sine of the angle, C is the cosine of the angle, i is the orbital inclination 
(97.58 degrees), e is the Earth spin axis tilt with respect to the ecliptic (23.442 degrees), u 


is the right ascension of the sun measured in the ecliptic plane and defined as 


u=u0+udot* time (2) 








with udot=0.985648 degrees/mean solar day; and w is the right ascension of the line of 


nodes defined as 
w=w0+wdot* time (3) 
where wdot is the nodal regression rate given by 
wdot=[[-9.9639(R/(R+h))*3/5]Ci]/(1-Ecc%2)“2 (degrees/mean solar day) (4) 


with Ecc being the eccentricity (zero for a circular orbit), h being the altitude of the 
satellite in orbit, 1 being the orbit inclination, and R being the radius of the Earth [Ref. 5]. 
The values for u0 and w0 are selected to vary Beta angle to test attitude dynamics 
assumptions for simulation runs. The reason for varying Beta angle is to generate 
different Earth to sun rotation angles which will be derived later. These rotation angles 


affect gravity gradient and aerodynamic disturbance torque inputs for the simulation. 
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Figure 1. Beta angle spatial relationships with Earth to Sun rotation angles 
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Figure 3. Beta angle spatial definition [Ref. 5] 


The key to the topic of a beta angle is the fact that the angles @,,, 0.,, and y,, in 
Figure 1 are dependent on beta angle and that these angles are vital to the coordinate 


transformation rotations needed to reference the Earth-oriented orbital frame to the 








Sun-oriented reference frame as discussed below. Using Figures 1 and 4 along with 


Napier's Circle relationships, the Earth to sun angles can be calculated as 


des = arcsin[cos B COS Ves | (5) 
Ges = arctan[ TS] (6) 
Wes = arcsin[cos B sin O¢5] (7) 


where f is defined as above and 
Ves = Vs + ol (8) 


where v, 1s a spatial factor necessary to orient the subsolar point in Figure 1 to the true 
anomaly v of the satellite. This is better seen by saying that the satellite starts its orbit at 


the equator (v=0) and proceeds around the Earth with angular rate w,. This means that 


The quantity v,, is also increasing at a rate equal to ,, but it starts (v., = 0) at the 
subsolar point instead of the equator; thus the need for the offset v,. Since the value of v, 
is dependent on the position of the subsolar point, it must be a function of the Earth's 
position in its orbit about the sun. This relationship can be developed using Figure 3 and 
some spatial relationship variables p, q, and t. A spherical triangle showing a portion of 
Figure 3 is given in Figure 5. Napier's circle relationships are used to solve for the angle 


q as 


g = arccos[—cos(180 — inc)cos e + sin(1 80 — inc)sine cos w] (10) 


10 





Now using spherical triangle relationships, p can be solved for as 


sin(180-inc)sin ud = (1 1) 


p =arcsin[ oe 


Using Equations (10) and (11) and the larger triangle of Figure 5, t can be solved for as 


t= arcsin[ SEG *y (12) 


Now using the smaller triangle and Napier's circle relationships of Figure 5 and Equation 


(12), v, can be solved for as 
v, = ¢—arcsin[tan B tan(90 — q)] (13) 


This offset becomes important later, because the Magnetic field model of the Earth is 


affected if this offset is not utilized. 





Figure 4. Napier's circle relationships for Earth to Sun rotation angle definition 
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Figure 5. v, spatial relationship generation 


C. REFERENCE FRAMES 

1. Earth-oriented Orbital Reference Frame 

This reference frame is the standard orbital coordinates reference frame that is 
used in most attitude dynamics publications. The positive X axis of the frame is oriented 
along the satellite's velocity vector with the positive Z axis located in the orbital plane in 
a nadir pointing direction (i.e. toward the center of the Earth). The positive Y axis is 
formed by taking the cross product of the Z and X axes. The Y axis is perpendicular to 
the orbital plane and along the negative orbit normal. 

2. Sun-oriented Orbital Reference Frame 

This reference frame is specially designed for this orbit application. In this 
reference frame, the negative Y axis is oriented in the direction of the sun and the plane 
formed by the Y and Z axes cuts through the center of the Earth. To orient the spacecraft 
in this reference frame, two axis rotations are performed from the Earth-oriented orbital 
reference frame using the Earth to Sun rotation angles discussed above. The first rotation 
is performed about the 3-axis or Z, axis of Figure 1 using the angle y,,. This forms a new 
set of axes called x1, x2, and x3. It can be shown that the x1-z1 plane and the y1-zl 
plane still cut through the center of the Earth since the z1 axis is still oriented toward the 
center of the Earth. Now a rotation about the x1 axis is performed using the angle 9,,. 
This rotation serves to now orient the negative Y axis in the direction of the Sun as in 


Figure 1. Now, since the x1 axis is perpendicular to the yl-z1 plane by definition, and 
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since the yl-z1 plane cuts through the center of the Earth before the rotation; it must still 
be cutting through the center of the Earth, which is the desired configuration. 

3. Body-fixed Reference Frame 

This reference frame is obtained by performing a standard 3-1-2 (Yaw-Roll-Pitch) 
Euler axis rotation in the Sun-oriented orbital reference frame. The rotations about all 
three axes will be small angles, so small angle approximations can be used in the 
definition of the equations of motion. 
D. SPACECRAFT BUS DESCRIPTION 

The spacecraft bus is shown in Figure 6. The main body of the bus has 
dimensions in inches of 19 x 19 x 30. There are three identical solar panels with the 
center panel mounted directly to the bus with the other two panels hinged to either side. 
All sun-monitoring equipment will be mounted on the top of the satellite. Figure 6 shows 
a rectangular box which will act as the equipment in this case for disturbance torque 
analysis. The momentum wheel and co-mounted Earth sensor are located at the back of 
the satellite toward the bottom. The Earth sensor assembly protrudes out the back of the 
satellite slightly. 

Inertia and center of mass properties for the satellite were extracted from 
Reference 4 and are listed below. The reference point (0,0,0) for all calculations is 
located at the center of the base of the satellite (i.e. on the center of the positive Z face at 


the bottom of the satellite). 


I. = 8.363 kg-m’ 
IL, = 9.284 kg-m? 
L, = 4.012 kg-m (14) 
I,y =-0.111 kg-m’ 
I, = 0.013 kg-m? 
I, = 0.048 kg-m’ 
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COM, = .000254 m 
COM, = -.005842 m 


COM, 


(15) 


453136 m 
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Figure 6. Spacecraft Bus [Ref. 4]. 


E. ATTITUDE DETERMINATION AND CONTROL COMPONENTS 


Table 1 gives a listing of control system components selected for the design. 


below. 


Each component will be explained thoroughly 


1. Sensors 


a. Earth Sensor 


the spacecraft attitude with respect to the Earth is of great 


Knowing 


importance in order to maintain constant communications with ground stations. Since the 


Earth appears to be rather large for near-Earth orbiting satellites, most sensors are 


designed to detect the Earth's horizon. 


14 









component size mass (kg) power (watts) comments 


Momentum 9.75 N-m-sec 8 3000 RPM 
wheel | +/- 10% 
Magnetic 15 A-m? 0.45 (x 3) 1.43 (x 3) dual coil 
Torque Rods 
Magnetometers | +/- 0.6 Gauss 0.25 (x 2) 0.04 (x 2) 
9 
8 
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5. 5 
Earth Sensor 1. 3.3 on momentum 
wheel initially 
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11.35 13.2 


The smallest amount of horizon sensing ambiguity is achieved by using 





Table 1. ADCS components 


the spectral region in the infrared near the 151m wavelength. Infrared detectors are not 
affected by night or the presence of a terminator. Furthermore, they are less susceptible 
to sunlight reflections off the surface of the spacecraft. However, visible light sensors 
have some advantages over infrared sensors including: lower cost, faster response time 
and higher signal to noise ratios due to increased visible radiated intensity. [Ref. 6] 

The main purpose of the Earth sensor for this satellite is to provide pitch 
error data which will be used to control pitch angle and keep the communications antenna 
(located on the Earth-facing side) nominally Earth-pointed. This will be accomplished by 
keeping the pitch/yaw plane in such an orientation that it cuts a great circle arc on the 
Earth. This guarantees Earth sensor data to the maximum extent possible. Due to the 
orbit configuration and attitude control approach selected, Earth sensor data will not 
always be available. Using Figure 7, one can see the relationships between the Earth 
sensor scan line, which is canted at a 45 degree angle up from the local vertical, and the 


surface edge of the Earth. It is easy to determine that if the Earth to Sun roll rotation 
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angle (@,,) is greater than 22 degrees in the negative direction, then the scan line of the 
sensor will not hit the Earth and sensing data will be unavailable. Current simulation 
results show that Earth sensor data will be unavailable for an increasing percentage of the 
orbit as Beta angle decreases. When Earth sensor data is unavailable, the magnetometers 
and an orbital elements model with a magnetic field model will be used to detect coarse 
pitch errors (magnetometer sensing error is approximately 5 degrees). 

The Earth sensor is mounted on the momentum wheel and protrudes out 
the positive pitch axis surface of the spacecraft (anti-sun face). Consideration was given 
to mounting the Earth sensor on the sun-facing side of the spacecraft for increased 
sensing data during eclipse when roll/yaw error signals are unavailable from the sun 
sensors, but solar array development considerations and the cost of an Earth sensor 
mounted separately from the momentum wheel can make this option unattractive. The 
addition of this extra Earth sensor may become necessary to provide sensing data in 


eclipse for low Beta angle applications, which will be explored in the simulations. 


direction of negative 


roll angle ee 
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scan line 
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Figure 7. Earth sensor coverage relationships 
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b. Sun Sensors 

Sun sensors have been used on nearly all spacecraft which have flown. 
Since the distance to the sun is rather large for Earth orbiting satellites, the angular radius 
of the sun is nearly independent of the type of orbit and for most applications, the sun can 
be treated as a point-source. Furthermore, the sun is bright enough to allow the use of 
simple equipment to discriminate among radiation sources with small power 
requirements. Generally, sun sensors are used for detecting sun presence, directing the 
solar arrays toward the sun and providing the relative sun vector for protection of 
sensitive equipment such as star trackers and other sensors. Sun sensors can have very 
narrow to very wide fields of view and resolutions from several degrees to less than one 
arc-second. The three basic types of sun sensors include: Analog sensors which provide 
a continuous output which is a function of sun angle; sun presence sensors which have a 
constant output when the sun is present; and digital sensors which provide a discrete 
output as a function of sun angle [Ref. 6]. 

High accuracy sun sensors will be used on the sun-facing side of the 
spacecraft to provide sun line roll/yaw error signals to the attitude control system for 
magnetic torque rod control. Coarse sun sensors will also be used on each face of the 
spacecraft in case a problem develops and the spacecraft's attitude with respect to the sun 
needs to be known for solar array pointing. These coarse sun sensors will also be needed 
to provide pointing control while reorienting the spacecraft from its Earth-oriented 
configuration during checkout to the sun pointing mode. 

c. Magnetometers 

Magnetometers are very popular attitude determination sensors and have 
been used widely on spacecraft. They are vector sensors because they provide both 
direction and magnitude of the ambient magnetic field. Furthermore, they are reliable, 
lightweight, and have low power requirements; they have no moving parts; and they can 
operate over a wide temperature range. Unfortunately, magnetometers are not very 


accurate inertial attitude sensors mainly due to the lack of an accurate model of the 
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Earth's magnetic field. Also, the Earth's magnetic field strength decreases as a function 
of the inverse cube of the satellite distance from the center of the Earth. This could result 
in residual spacecraft magnetic effects dominating the field strength, generally limiting 
the use of magnetometers to spacecraft below 1000 kilometers (note: missions have 
successfully used magnetometers above this altitude) [Ref. 6]. 

The magnetometer consists of a magnetic sensor and an electronics unit 
for transforming sensor measurements into a usable format. There are two main types of 
magnetometers: quantum magnetometers, which utilize atomic properties; and induction 
magnetometers, which are based on Faraday's Law of Magnetic Inductance. The latter 
group is more popular and will be discussed further here. [Ref. 6] 

Faraday's law states that an electromotive force (EMF) is induced in a 
conducting coil when it is placed in a time-varying magnetic flux. There are two types of 
induction magnetometers called search-coil and fluxgate magnetometers. Search-coil 
magnetometers contain a solenoid coil of N turns surrounded by a ferromagnetic core 
with magnetic permeability 1, and cross-sectional area A. The EMF induced in the coil 
when placed in a magnetic field produces a voltage as a function of the time varying 
component of the magnetic field which is parallel to the solenoid. This type of 
magnetometer is mainly used on spin-stabilized spacecraft to provide precise phase 
information. The fluxgate magnetometer uses two coils and two saturable cores to detect 
the presence of any ambient magnetic field by measuring harmonic frequencies formed 
by the time-varying field. 

A pair of two-axis magnetometers will be used to detect the ambient 
magnetic field the spacecraft is encountering. The data will be used primarily to calculate 
electrical currents needed to generate desired control torques with the magnetic torque 
rods and also in a backup roll to detect course attitude errors as discussed in the Earth 
sensor section above. The magnetometers will be mounted in such a way that each 


element detects a component of the magnetic field along each major body axis. These 
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signals will then be averaged to determine the magnetic field strength and direction. This 
configuration will provide some backup in case of failure in one of the components. 
2. Actuators 

a. Momentum Wheel 

Angular momentum storage devices are used on spacecraft for several 
purposes including: the addition of stability to counter disturbance torques, to provide a 
variable momentum to allow operation at one revolution per orbit for Earth oriented 
missions, to absorb cyclic torques, and to transfer momentum to the spacecraft for 
slewing maneuvers. The principle of momentum for a spinning wheel is derived from the 


relationship 
h=/o (16) 


where J is the moment of inertia of the wheel and wis the rotational velocity. 
Momentum wheels are designed to provide a variable momentum ‘bias’ or 


nonzero momentum which can be used for storage capability about the rotational axis of 
the wheel which is generally fixed with respect to the spacecraft. The same momentum 
can be achieved with a large inertia, slowly rotating wheel and a small inertia, rapidly 
rotating wheel. Design tradeoffs usually favor the smaller wheel for the weight savings, 
but the increased wheel speed can lead to greater wear on the bearings which may shorten 
the wheel lifetime. 

Torque motors using either direct or alternating current (DC or AC) are 
used to transfer momentum between the wheel and the spacecraft body. The DC motors 
are generally more efficient and provide high torques at low speed, but they do require 
the use of brushes which can wear out. AC motors do not have this drawback, but they 
do require gearing systems due to their high speed. Also, they are less efficient and 


provide a lower torque. Tachometers are used to measure wheel speed for attitude control 


purposes. 
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The main purpose of the momentum wheel on this spacecraft is to provide 
gyroscopic stiffness along the pitch axis once the spacecraft is operational. It also serves 
to control pitch angle by changing wheel speed to generate the necessary control torques. 
The amount of time that the momentum bias vector needs to keep pointing accuracy is 
assumed to be 35 minutes; the maximum eclipse time. Momentum wheel size is listed in 
Table 2 for several pointing accuracies over this maximum eclipse period assuming worst 


case disturbance torques are present throughout the eclipse. This assumption provides for 


some built in safety margin for the sizing. 


accuracy after eclipse period (degrees) momentum wheel size (N-m-sec) 





Table 2. Momentum wheel size versus gyroscopic stiffness pointing accuracy 


The amount of torque output required is driven by Iyy and pitch angular 
acceleration (6), which is a function of Beta angle. Equations for 6 are developed in 
Reference 4 and are listed here as 


w2K sink 
8 max = (1-0 cos)? (17) 
where 
4cos*PsinB 
nt cache ao l 
(14sin?B)? Ge) 


20 








_ 1-sin7pB 
l+sin7B 


(19) 


148602 -1 


€ = arccos[——>——] | (20) 


The maximum momentum change required from the wheel is given by 


. 2 
AH yneet = Omax — Ominhy =(Aq-Sin B\Oohy =FFoohy (21) 


The maximum torque required from the wheel is then determined from 


T max = Lyy8 max (22) 


A wheel with a momentum of 9.75 N-m-sec at 3000 rpm with a torque rating of 0.025 
N-m was selected. This provides a pointing deviation slightly less than 0.5 degrees over 
the maximum eclipse period assuming worst case disturbance torques are acting on the 
spacecraft throughout this period. Initial simulation results indicate that the maximum 
pitch angular acceleration can be shown to be approximately 2 x 10° rad/sec*. This 
means there can be an extremely large Iyy before wheel torque becomes a concern. 

b. Magnetic Torque Rods 

Magnetic coils or electromagnets are used to generate torques necessary to 
control spacecraft attitude and momentum. This could include use for overcoming 
disturbance torques and compensating for residual magnetic effects. The magnetic 


moment generated by a loop of wire possessing N turns is given by 
m = N/An (23) 


where N is the number of turns, I is the current through the loop, A is the area enclosed 


by the loop and n is a unit vector normal to the plane of the loop. The positive sense of 
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the direction of the magnetic moment is determined by the right hand rule. The fingers of 
the right hand are curled around the turns of wire in the direction of the current and then 
the right thumb points in the direction of the magnetic moment. 

The magnetic dipole moment depends on what type of material is enclosed 


within the turns of wire and is expressed by 
M=pm (24) 


where 1 is the permeability of the core material. The permeability of free space 1, has 
the value 4 x 10° N/A’. Thus, it is clear that the values of core permeability, coil 
configuration, N and A, and the current values must be selected appropriately for the 
given application. Magnets on most satellites have used an "air" core. To generate 


magnetic control torques the following relationship is used 


where M is the magnetic dipole moment and B is the magnetic field expressed in body 
coordinates. | 

The magnetic torque rods will be used for momentum wheel desaturation 
and roll/yaw repositioning to zero out sun sensor error signals. The change in momentum 
required to reorient the momentum vector by 1 degree per day is approximately 
0.122 N-m-s. This relates to a constant torque of 1.41 x 10° N-m throughout the day (or 
2.82 x 10° N-m on a 50% duty cycle) . This value is much smaller than the maximum 
disturbance torque, so the latter will be used to size the magnetic foraue rods. Assuming 
a worst case magnetic field strength of 1.5 x 10” tesla, the required torque rod moment is 
approximately 2.45 A-m’. Magnetometer sensing data will be used to provide control 


signals for electrical current generation which control the direction and amount of torque 
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output. The selected torque rods have two coils, which provides redundancy in the 
system, and have a magnetic moment of 15 A-m’. 

c. Thrusters 

Since no orbit control is considered to be necessary, a thruster system 
was omitted due to a great increase in weight, complexity and cost. 
F. SYSTEM PERFORMANCE 

Once the injection vehicle has attained the proper orbit, the momentum wheel will 
be spun up with the injection vehicle providing torque to counter the increase in 
momentum. It will take approximately 7 minutes to spin up the wheel if injection vehicle 
control is available. In the event the launch vehicle cannot provide stabilization, the 
wheel will be spun up in accordance with the capture control section in Chapter IV. 
When the wheel has achieved speed, the separation system will activate (assuming 
injection vehicle control) and separate the satellite from the injection vehicle. Once clear, 
the solar arrays will deploy and operations will begin to establish a stable nadir-pointing 
attitude. Worst case tip off rate in pitch is expected to be approximately 5.5 degrees per 
second and worst case nutation angle (due to roll tip-off torque) is expected to be 5 
degrees. The momentum wheel is designed with sufficient momentum storage to absorb 
the former and the magnetic torque rods will be used to damp out the latter. 

When the satellite is stable in a nadir-pointing configuration, system checks will 
be performed. Next, ground commands will be used to initiate the reorientation of the 
momentum bias vector toward the sun with torque rod control. Following this, the 
attitude control system will maintain pointing accuracy of 0.5 degree in roll and yaw and 
5.0 degrees in pitch. 

The allowable speed region for the wheel will be 3000+/-300 RPM. Wheel 
desaturation will be performed as needed using the magnetic torque rods. Initial 
simulation results show that wheel speed takes several orbits to reach desaturation limits 
and that the magnetic control law is easily capable of accomplishing the desaturation. 


Other simulation results from MATLAB/SIMULINK also indicate that the momentum 
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bias vector is capable of limiting roll/yaw drift to .2 degrees throughout the maximum 
eclipse period (35 minutes) with no control inputs. If the proper alignment is achieved 
before eclipse entry, then the spacecraft should be oriented at visible sunrise so that the 
sun is within the field of view of the sun-monitoring instruments. 
G. SIMULATION 

An old SIMULINK program, developed by Professor Leonard at the Naval Post 
Graduate School, Monterey, California has been upgraded and tailored for this particular 
application. It includes simulation blocks for each disturbance torque (solar pressure, 
gravity gradient, and aerodynamic) and a simple simulation for the Earth sensor and sun 
sensors. From these inputs, spacecraft torques will be calculated and the spacecraft's 
dynamic properties will be simulated. A separate MATLAB program will be used to 
control the SIMULINK program and create output displays. SIMULINK and MATLAB 


block diagrams and control codes are listed in the appendix. 
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lI. SPACE ENVIRONMENT AND ITS AFFECT ON THE 
SPACECRAFT 


A. EARTH'S MAGNETIC FIELD DEVELOPMENT 

The derivation of the Earth's magnetic field is very complex for both magnitude 
and direction. However, the field can be approximated by a large, tilted dipole located at 
the Earth's center with the ‘south’ end of the dipole located in the northern hemisphere. 
The location of the magnetic poles is constantly drifting at approximately .014 degrees 
per year [Ref. 6]. More accurate models of the field are available but they require vast 
computing resources to simulate; and they are not required for use in this feedback 
control application. Thus, the dipole approximation will be used. 

The plane which is perpendicular to the dipole is called the magnetic equator and 
this is where the field is weakest in magnitude. Figure 8 shows how the madenetic field 
strength varies with altitude at the magnetic equator. It can be shown that the strength of 
the dipole field decreases as a function of the inverse cube of the distance from the dipole 


center (i.e. the center of the Earth) [Ref. 6]. 
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Figure 8. Magnetic Field Intensity at Equator in (nano-Tesla) [Ref. 6] 
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As magnetic latitude increases from zero to 90 degrees, the field strength at a given 


altitude doubles as shown in Figure 9. 
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Figure 9. Relative Intensity of Magnetic Field as a Function of Latitude [Ref. 6] 


For magnetic control, the Earth's magnetic field needs to be expressed in 
coordinates relative to the three primary body axes of the spacecraft. The magnetic field 


components expressed in the Earth oriented orbital reference frame are: [Ref. 7] 


Byo = o31CeSiCo. — SaSpSe — CiCaCuSe] (26) 
Byo = a [-CeCi — CuSiSe] (27) 
Bro = a3 L2CeSiSo + 2CaSpSe — 2CiSaCuSe]. (28) 
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Where C is the cosine of the angle, S is the sine of the angle, i is the orbit inclination, R is 
the radius of the satellite orbit, Mis the Earth's Dipole strength (7.943 x 10’ Nm’/Amp), 
¢ is the dipole tilt with respect to the Earth's spin axis (11.7 degrees), 1 is the Earth's 
rotation with respect to the ascending node (varies 360 degrees in one day), and o is the 
orbit angle with respect to the ascending node (i.e. the angular position of the satellite 
measured counterclockwise around the orbit from the ascending node). 

To express these magnetic field components (B,,, B,,, and B,,) in the spacecraft 
body axis reference frame, two coordinate transformations are needed. The first 
transforms the magnetic field from the Earth-oriented orbital reference frame to the 
sun-oriented orbital reference frame, while the second is a standard 3-1-2 Euler rotation 
which leads to the body-fixed reference frame orientation. Once the field components are 
represented in body-fixed coordinates, they can be used in the equations of motion for 
both disturbance torque analysis and magnetic control torque inputs. 

B. TORQUES WHICH AFFECT THE SPACECRAFT 

To be able to analyze the motion of the spacecraft on orbit one must be able to 
model all torques which act upon the spacecraft whether they are generated by the 
environment or by the spacecraft attitude control system. These torques can either be 
categorized as cyclic, varying sinusoidally over an orbit, or secular, accumulating with 
time and not averaged over an orbit. Numerical integration of Euler's equations requires 
that all torques be modeled as a function of time and the spacecraft's attitude. For this 
analysis, the major environmental disturbance torques to be considered are: aerodynamic 
drag, solar pressure, residual magnetic dipole and Earth's gravitational. These 
disturbance torques would quickly reorient the spacecraft if not resisted in some way by 
the attitude control system. The active control torques to be used are generated by the 
momentum wheel and magnetic coils. 

1. Disturbance Torques 

Before sizing the ADCS components, one must know the maximum possible 


disturbance torques that the spacecraft might experience throughout its lifetime. As 
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discussed above, the major categories of disturbance torques are: aerodynamic, solar 
pressure, residual magnetic dipole and Earth's gravitational. The type of torque which is 
the largest for this type of low altitude application is the aerodynamic torque. Also, since 
there are several techniques available to limit the effects of residual magnetic dipole 
disturbances, this category will be ignored for the purpose of simplicity in the simulation. 
Using References 6, 8 and 9, each disturbance torque is analyzed to determine its 
magnitude and which axes of the spacecraft it affects. 

a. Aerodynamic 

For lower altitude orbits, the interaction of the spacecraft with the Earth's 
upper atmosphere causes a disturbance torque about the spacecraft center of mass. The 
force due to the impact of atmospheric molecules on the spacecraft surface can be 
modeled as an elastic impact without reflection. The incident particle's energy is 
generally completely absorbed [Ref. 6]. The incremental force dF on a surface element 


area dA with outward normal N is given by 
dF =-5CppV*(N-V)Vdd (29) 


where V is the unit vector along the translational velocity of the individual surface 


element relative to the incident stream of particles, Cp, 1s the coefficient of drag (assumed 
to be 2.0 if no measured value is available) and p is the atmospheric density. To 


determine the total aerodynamic torque, the individual forces must be integrated as 


follows 


T aero = [Rer-ce x dF (30) 
where R,,., 1s the vector measured from the center of mass of the spacecraft to the center 


of the individual incremental areas (i.e. the distance to the center of pressure for that 


individual area). 
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The aerodynamic disturbance torque is a fixed quantity for spacecraft 
which are oriented toward the Earth, but varies for those which are inertially referenced 
such as this spacecraft. The velocity vector for the spacecraft can be expressed in 


Earth-oriented coordinates as 
V sc = VXo (3 1) 


To take the necessary cross products for determining the aerodynamic torque, this 
velocity vector must be expressed in body coordinates by performing the two required 
coordinate transformations using the Earth to sun angles and the sun to body Euler 
angles. Since the vector r,,., for each individual piece of the satellite is known in body 
coordinates, it is easy to simulate the cross products and thus, the aerodynamic torques 
for each individual area of the spacecraft. These individual torques are then added to 
determine the overall aerodynamic torque. 

From Reference 9, the worst case atmospheric density is approximately 
2.0x 10" kg/m*. Using a coefficient of drag of 2.0, orbital velocity of 7585.1 m/sec and 
the effective surface area for each individual satellite panel the aerodynamic force on 
each face of the satellite can be calculated and the individual cross products can be 
determined. The individual, time varying computations are performed in the SSMULINK 
program and will not be listed here. 

b. Solar Pressure 

Radiation incident upon the body of the spacecraft produces a small 
disturbance torque about the center of mass of the spacecraft. The surface of the body is 
subjected to a radiation pressure (force per unit area) which is equal to the vector 
difference between the incident and reflected solar radiation energy. Since the intensity 
of the solar radiation is a function of distance from the sun, spacecraft altitude with 
respect to the Earth does not affect solar radiation pressure for Earth-orbiting spacecraft. 


The mean momentum flux acting on a surface normal is determined from 
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p= (32) 


where F. is the solar constant (1358 W/m’) and c is the speed of light (3.0 x 10° m/sec). 
This flux value is constant except for a small periodic variation due to the eccentricity of 
the Earth's orbit about the sun [Ref. 6]. 

For most applications, the solar force can be modeled adequately by 
assuming the incident radiation 1s either absorbed, reflected specularly, reflected diffusely, 
or some combination of the above. Specular reflection means that the incident vector is 
reflected at the same angle as it came in to the surface, while diffuse reflection means that 
the incident vector is reflected in all directions. The differential force imparted by the 


particles which are totally absorbed 1s given by 
d¥ obs = —PCacos9SdA (33) 


where C, is the absorption coefficient, S is the unit vector from the spacecraft to the sun 
and 0 is the angle between S and N, the unit vector normal to the differential surface. This 
last angle must be less than 90 degrees for a pressure force to exist on that area. The 


differential force for the specularly reflected radiation is given by 


where C. is the coefficient of specular reflection. For a surface which reflects diffusely, all 
reflection angles must be considered and integrated to form the total differential force for 


diffuse reflection as 
dF gig = PCA(-4 cos8S — cos8S)dA (35) 


where C, is the coefficient of diffuse reflection. If all types of reflection play a part, then 


the total differential force is given by 
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AF total = —P (ta —Cs)S + 2(Cscos8 + 5Ca)N]cos 6dA (36) 


where C,+ C,+ C, = 1 [Ref. 6]. 
The solar torque is then determined by taking the integral of all cross 
products between each differential solar force and its individual distance offset from the 


center of mass of the spacecraft [Ref. 6]. 
T solar = | Mep-cg X AF total (37) 


Since it is difficult to perform this integral, it is not unusual to break down the spacecraft 
to simple geometric elements and then estimate each of the individual torques for each 
element. This spacecraft will be treated in a similar manner. Since only the solar panel 
side of the spacecraft will be pointed at the sun during normal operations, we must 
consider that side only for solar pressure (i.e. the angle between S and N for the other 
sides is 90 degrees or the face is not illuminated). This means no solar torques about the 
Y axis will be generated. The surfaces on this side of the spacecraft include the solar 
panels and the sun-facing side of any monitoring instruments [Ref. 6]. 

For the solar arrays, the total surface area is 1.094 m’ and it is assumed 
that C,= 0.8, C,=0.1, and C,=0.1. The center of pressure for the array assembly 1s 
located at (x, y, z) = (0.0, -9.5, -14.9) inches or (0.0, -0.2413, -0.3785) meters using the 
body coordinate reference listed earlier. Utilizing Equations (15) and (37), this will cause 


a torque about the negative X and Z axes equal to 


Tspx = —4.31x107’ N-m (38) 
T'spz = —1.47x107? N-m 


The sun-monitoring equipment is assumed to have a surface area of 0.061 


m? and it is assumed that C,= 0.05, C,= 0.8, and C, = 0.15; with a center of pressure 
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located at (x, y, Z) = (0.0, -9.5, -32.25) inches or (0.0, -0.2413, -0.8192) meters. This will 


cause a torque about the positive X axis and the negative Z axis equal to 


T egpmntz = —1.33x107!9 N- m 


This yields a total solar torque for each axis of 


T solarx = —2-39x1077 N—m (40) 
T solarz = —1.60x10-? N-—m 


c. Gravity Gradient 

Any non-symmetrical body in orbit will experience a gravitational torque 
because of the variation of the Earth's gravitational force over the object. This is due to 
the non-uniformity of the inverse square gravitational field. To develop the relationships 
for the gravity gradient torques it is assumed that the moments of inertia are known and 
that the spacecraft is orbiting a non-realistic, spherical Earth. Using Reference 8, the 
differential force acting on the spacecraft mass element dm located a distance R from the 


center of the Earth is given by 


eRdm 
dF = —"" = (41) 


where pe 1s the Earth's gravitational constant and R is the distance from the center of the 
Earth to the individual incremental masses. Using Figure 10, R can be expressed as 


R=R,-r (42) 


where R, is the distance from the center of the Earth to the center of mass of the 


spacecraft and r is the distance between the spacecraft center of mass and the individual 
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incremental masses. Substituting this relationship into the above equation and using the 


fact that r/R, << 1 to expand R” in a Taylor series expansion yields 


e(Ro-r)di ‘Ro 
dF x a? + 3 (43) 


The torque about the center of mass is obtained by integrating each of the — 


individual incremental forces and noting that 
fr dm = 0 (44) 


because r is a vector measured from the center of mass of the spacecraft. 





Figure 10. Gravity gradient spatial relationships [Ref. 8] 


The final result is 





re le xR,)(r-Ro)dm (45) 


3 
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The vector R, , which is expressed in the Earth-oriented reference frame as R, K, can be 
written in body coordinates by performing a coordinate transformation from the 


Earth-oriented to the body coordinate system. The result is 


Ro = Ro Ai body + Bj body + Ck body) (46) 
where the expressions for A, B and C are given by 

A = Sdes[COSy + SOSOCY] — Co es[SOCO] 

B= Shes[CyCo] + Coes[SO] (47) 

C = Sdes[SOSw — SOCOCW] + Coes[COCO] 


where C is the cosine of the angle and S is the sine of the angle. Substituting Equation 
(46) into Equation (45) and integrating yields the gravity gradient torques about each axis 
in body coordinates (ignoring cross-product terms) as 
| a 
T gox = = Uz — Iyy)BC 
3ue 
3ble 
P'ggz = Ty Ux — by )AB 
d. Residual Magnetic Torque 
Through the use of non-ferrous materials and proper wiring 


configurations, the residual magnetic dipole of the spacecraft can be greatly reduced and, 


in most cases, neglected. There can be no specific determination of the residual magnetic 
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dipole until the spacecraft is actually built and tested. For simulation purposes, residual 
magnetic torques will be ignored in this application. 
e. Total Torque 


Initial simulation results show that the estimated worst case torques for 


each axis are expected to be approximately: 


roll 2.0x 10° N-m 
pitch 9.0 x 10° N-m (49) 
yaw 2.0 x 10° N-m 


Since the yaw and roll torques are expected to be the largest, they will be used to size the 


ACS equipment. 
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IV. EQUATIONS OF MOTION AND ATTITUDE CONTROL LAWS 
A. EQUATIONS OF MOTION 


From Reference 8, the basic moment or torque equation is 


To. =H) + © x He (50) 
where H, is given for a system with a momentum storage wheel by the equation 


H- = 10 + Dyheet (51) 


where I is the inertia matrix of the satellite and @ contains the body rotational rates of the 


satellite given by 
x = COb — CHSOW —  on(COSy + SHSOCy) 
Oy = 6 — Shy — © onCoCy (52) 
wz = SOb + CHCOY — © on(SOSy — ShCOCw) 


where 6, 6 and y are the angles used in the 3-1-2 Euler angle rotation from the 
sun-oriented orbital reference frame to the body reference frame. The term @,, is the 
angular velocity of the sun-oriented reference frame. Assuming the time rate of change 


of Beta to be zero, it is given by 


cos700 9 (53) 


on = Ves ® Fa 2s 
cos“vsinB 


Equation (50) can be written in components as 
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Ty = H,+@yH;-o-Hy 
7. = H, + oxHy — oH, 


A special case to consider occurs when the rotating axes x, y and z coincide with the 


principle inertia axes of the spacecraft. Then Equation (51) simplifies to 
Hy = Ty x + Aapheetx 
Ay = Iy@y + hyheely (55) 
Az = [7202 + Myheelz 


These equations can be substituted back into Equation (54), noting that for this 


application there is only wheel momentum about the y axis, to form 
Ty = Ip, x + My@zUzz — Ly) — Ozh wheely 
Tol ld + Ops hele) eh ed (56) 
Tz = 17202 + OxOy(Lyy — Lex) + Oxhwheely 


These equations can be manipulated, using substitutions from above and noting that the 
wheel momentum is oriented along the negative y axis, to form the linearized moment 


equations for each major axis 
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Tox + Tay = Lex — © onW — ® on) + Cy — Izz)(@ on + (ond) + A(W + @ ond) 
Tey + Tay = Iyy(6— Gon) —h | (57) 
Tez + Taz = Tez( VU + O ond + © ond) — Cy - Tx)(@ ond — on) —h(b - © ony) 


where T, constitutes the control torques and T, constitutes the disturbance torques. 


B. CONTROL FLOW 
Figures 11 and 12 present the control flow for the spacecraft attitude control 


system and Earth sensor/Pitch control. 
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Figure 11. Spacecraft attitude control flow 
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Figure 12. Wheel/Pitch control flow 





C. ATTITUDE CONTROL LAWS 

1. Satellite Capture After Booster Separation 

Once the launch vehicle achieves the proper Earth-oriented configuration, the 
satellite will be separated. Two possible separation scenarios will be considered and 
simulated. In the first case, the momentum wheel will be spun up with the launch vehicle 
providing attitude stabilization until the wheel reaches optimum speed. Upon separation, 
only disturbance torques and torques generated by the separation mechanism will cause 
the satellite to tumble. For the second case, it is assumed that the launch vehicle cannot 
provide attitude stabilization. The satellite will be separated from the launch vehicle and 
the momentum wheel will be spun up at the maximum torque rating. This procedure, 
along with any torques generated by the separation system and disturbance torques, will 
cause the satellite to spin at an increasing angular rate about the pitch axis until the wheel 
is up to speed. It 1s assumed here that only one spring failure will occur among the four 
springs 1n the separation system. This single spring failure will cause the satellite to 
tumble if the wheel is not up to speed, and will cause nutation or pitch error if the wheel 


is up to speed. All tumbling motion needs to be removed so the satellite can be stabilized 


40 








in an Earth-oriented configuration for system check out and subsequent transition to the 
sun-oriented mode. 


Reference 2 suggests the use of the following magnetic control moments 


M, = _KB, 
My, =—KBy (58) 
M, = ~KB, 


where K is the control gain and the B; components are the time derivatives of the 
magnetic field relative to the body axes of the satellite. These quantities can be generated 
by taking the derivative of the magnetometer outputs [Ref. 2]. 

This theory functions on the premise that the kinetic energy of a body, due to the 
rotation about its center of mass, can be reduced; which eliminates all angular motion 
between the satellite and the magnetic field. This results in a terminal motion where the 
satellite rotates twice per orbit about the pitch axis (once per orbit relative to the Earth- 
oriented orbital reference frame) [Ref. 2]. 


The time derivative of the rotational kinetic energy is given by 


Beat (59) 


where t is the external torque acting on the body and o is the rotational velocity vector. 


A torque generated by the magnetic control torque rods is given by 


Now the derivative of the magnetic field can be estimated by (assuming @>>o,) 
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Be Bx@ (61) 
Combining the above equations, yields 
dkE =~ B-M (62) 
Implementing the magnetic moment equations 1n Equation (58) yields 


Oe = -KB? =-K|Bx |? (63) 
This quantity is clearly less than or equal to zero which is desired to remove rotational 
energy. Equation (63) shows that this law cannot reduce the kinetic energy associated 
with the component of w parallel to the magnetic field. Fortunately, the orbital motion of 
the satellite will ensure that the direction of the magnetic field relative to the body is 
continuously changing [Ref. 2]. 

Results are shown in Figures 13 through 15. A control gain of 1 x 10’ (N-m-sec’ 
/tesla’) was used to maintain torque rod saturation throughout the capture procedure. 
Figure 13, shows the results of a simulation where the wheel has already been spun up 
and the initial pointing error between the orbit normal and the spin axis is 20 degrees. 
This was done to show that the system can remove the initial nutation due to the inertial 
body rotational velocities and also remove the pointing error between the orbit normal 
and spin axis. Visible is the initial chaotic reaction to the separation where the rapidly 
changing initial response is due to nutation. Once all body rotational rates have been 
neutralized, the system begins to remove the pointing error. This particular stmulation 
took approximately 2 orbits to remove the pointing error between the orbit normal and 


the spin axis. 
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orbit normal/spin axis angular error in degrees vs time 
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Figure 13. Initial capture pointing error to orbit normal. Momentum wheel 
spun up before separation 


Figure 14, displays the simulation results for the case where the satellite 
momentum wheel has not been spun up and the satellite is ejected with initial rotational 
velocities about the x and y axes of approximately five degrees per second. As is clearly 
apparent, the satellite tumbles through several revolutions before the momentum bias 
overcomes the body angular momentum magnitude, which is being continuously 
decreased via the control law. This simulation shows that the satellite went through 30 
complete rotations before it settled in the proper configuration. There is, however, an 
inherent problem with this approach to capture when the wheel is not up to speed at 
separation. The spin axis will always align itself with the orbit normal, but there is the 
possibility that the satellite will be turned about the yaw axis by 180 degrees (momentum 
vector pointed in opposite direction to that desired). In this case, ground command will 
be necessary to precess the momentum vector through that 180 degree swing to get into 


the proper configuration. Figure 15 is an expanded view of Figure 14. 
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Figure 14. Initial capture pointing error to orbit normal. Momentum wheel 
spun up after separation. 


2. Momentum Wheel Control 
Wheel momentum will be controlled using the proportional/differential control 


law 
h=—Ko(t90 + 8) (64) 


where @ is the pitch error signal provided by the Earth sensor or magnetometers and 
Ko and t¢ are control gains [Ref. 8]. Substituting Equation (64) into Equation (57) yields 


Ignoring the wheel desaturation control input T, for now and considering the last term on 
the right hand side to be part of the disturbance torques (this term is on the order of 10° in 


magnitude or it can be canceled by feeding forward an estimate of the signal) yields 
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orbit normal/spin axis angular error in degrees vs time 
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Figure 15. Exploded view of Figure 14 showing final pointing angle. 
Taynew = Ly + Koto + Ko0 (66) 
Transforming to the Laplace domain 


Taynew(S)  s20(s)}+Koto/Iyys0(s)+Ko/Iyy0(s) 





This is the transfer function for a second order damped oscillator with natural frequency 


On= |= (68) 


g=F Ii (69) 
The design point is for a critically damped system where E=1 and then 
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I 
t= 2.) (70) 


If K, is selected as 0.36 kg-m’/sec’ and t, is selected as 10 seconds, the natural frequency 
is 0.19 radians/second with a settling time of approximately 40 seconds, assuming 
settling time is defined as four times the time constant Tp. 

System response to an initial pitch error of five degrees is shown in Figure 16. 
The system is reacting as a critically damped oscillator with a settling time of 


approximately 40 seconds, as designed. 
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Figure 16. Momentum wheel control response to initial pitch error. 
An interesting pitch response occurs when Earth sensor data is unavailable and 


only magnetometer pitch error is fed back for control. Since pitch rate is no longer 


available, Equation (67) becomes an undamped oscillator which will grow over time with 


46 








pitch disturbance torque input. Fortunately, Earth sensor data is unavailable for only a 
short period of time, so the pitch response does not grow too large as in Figure 17. If this 
oscillation becomes a problem for the overall system response, then some rough 


estimation of pitch rate will need to be generated to provide some damping. 


pitch response after magnetometer data feedback 
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Figure 17. Pitch response with magnetometer pitch angle feedback only 


3. Momentum Wheel Desaturation Control 

If a secular disturbance torque is experienced about the pitch axis, momentum 
wheel speed will grow without bound to counteract the disturbance torque [Ref. 8]. Once 
the momentum wheel's speed has deviated far enough away from its nominal speed (3000 
rpm), something must be done to return the wheel to the correct speed. That is the 
purpose of the momentum wheel desaturation control. The portions of the magnetic 


control moments for wheel desaturation are given by 


47 








M, = —K,,BzAh 
My =0 (71) 
M, = K,B,Ah 


where Ah is the overall change in the wheel's momentum from its nominal value and K, is 
the control gain for desaturation. After the cross product is taken with the magnetic field, 


the control torques are calculated to be 


‘- = -KyBxByAh 
Ty = K,(Bz +Bz)Ah (72) 
Tz = —K,ByBzAh 


Of note, this control law generates the required torque about the y axis for wheel 
desaturation, but it also creates disturbances about each of the other axes which must be 
accounted for. Also, the control gain, K, , must be a negative number. This can be seen — 
by the need for a negative y axis torque if Ah is a positive quantity and vice versa. The 
control gain K, is sized to provide torque rod saturation when Ah is ten percent of nominal 
wheel momentum. This leads to a control gain of -1.538 x 10°(sec" - tesla”). 

Momentum wheel desaturation control will be shut off when sun sensor data is 
unavailable because the control law creates disturbances for the x and z axes which cannot 
be dissipated while there is no sun sensor data. 

4. Control for Transition from Earth-oriented to Sun-oriented Attitude 

The control law for sun-oriented control was developed first with the hope that it 
could be used without modification to accomplish the transition from Earth-oriented/ 
checkout operations to sun-oriented operations. The orbit configuration will be adjusted 


in order to vary the Earth to sun angles o,. and w., and test the control law at several 
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values of Beta angle. It was not known whether the sun-oriented control law could 
handle such large angular errors since the control law was based on small angle 
assumptions). | 

Reorientation simulations were only performed for a 1:30 orbit and a 3:00 orbit. 
It was assumed that if the reorientation worked for these initial conditions, then it would 
work for the 4:30 and 6:00 orbit because the transition angles will be smaller and 
disturbance torques will be smaller due to the lower Earth to sun angular values. 

Figure 18 displays the results of a reorientation simulation from an initial sun 
pointing error of approximately 72 degrees. This constitutes a reorientation for a 1:30 
am/pm orbit with average Beta angle of 18 degrees. It takes several orbits to perform the 
reorientation as expected. Figure 19 shows a reorientation simulation for a 3:00 am/pm 
orbit. The raised dotted lines at the bottom of the plots represent the time periods when 
the satellite is not in eclipse. If the value of the dotted line is one, then roll/yaw attitude 
sensing data is available. The eclipse control also serves to cut off momentum wheel 
dumping control when the satellite is in eclipse. 

5. Sun-oriented Attitude Control 

The pitch axis control torque will be generated as discussed in the momentum 


wheel and wheel desaturation control sections. The magnetic moment equations are 


given as 
M,, = —K,BzAh 
My = K, Bx + K2(Bzo — Bp) + K3Bzy (73) 
M, = K,B,Ah 


where K, and K, are used for precession control, while K, is used to control nutation. 


Reference 10 used the control law 
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For small angles, since no measurement for yaw existed; By can be approximated by 
B, 
By = (Bzo —- Bx) - Oo(—— +20Bx) (75) 


In the orbital reference frame, B, is constant for an Earth-pointing satellite. For a 
sun-pointing satellite, it is not constant, so the control law from Reference 10 should not 
be used. Instead, the y axis torque rod control moment 1s selected as the second equation 
of Equation (73), which is quite stmilar to Reference 10. The additional K,B,w term has 
been added because there is a yaw error measurement which is not available when using a 


horizon sensor, which is the standard sensor for Earth-pointing satellites. 


pointing error in degrees vs time in minutes; reorientation contro! 
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Figure 18. Earth to sun reorientation (1:30 orbit) 
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pointing error in degrees vs time; sun-oriented control 
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Figure 19. Earth to sun reorientation (3:00 orbit) 


When the M x B cross product is performed the following control torques are generated 
T, = K,BxBzb + K7(B2b — ByBzW) + K3B2y — KyBxByAh 
Ty = Ky(Bz + Bz)Ah | (76) 
T, =—K Bz — K2(BxBz — Bz) — K3B,Bzy — KyByBzAh 


When these control torques are substituted into Equation (57) above, the total torque 


equations become 


la Txx — — K>B2$- [Cx — Ly + Izz)® on +h- — K>B,Bz|w + KyByB,Ah 
+[(hyy — Lec) on + h® on — K,BxBz]o - nx® on + K3B2]y 


T ay = ly + Kot90 + Ko® — Ku(Be + Bz)Ah— I on (77) 
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Tz = Tez + [Ulex — Ty + Le.) on —h + Kp BxBz\b — KyB2yy + Ky,ByBzAh 


Figures 20 and 21 show the reason for selecting the control laws as listed. In Figure 20, a 
positive ¢ error is shown which requires a change in momentum in the positive z 
direction. This can be generated by the (-K,B,’ $) portion of the z axis magnetic torque in 
Equation (76). If K, is selected as a negative number, this generates a positive torque as 
required to achieve the desired change in momentum. Figure 21 shows what happens 
when there is a negative y error. This error requires a change in momentum in the 
positive x direction which is generated by the (K,B,” w) portion of the x axis control law 
in Equation (76). Of note, the portions of the x and z axes control laws which contain 
B,B, create disturbance torques about their respective axis. However, over the course of 
one orbit, these disturbances should average out to zero. This averaging could be 
precluded by the effects of eclipse which would cut off attitude data for a portion of the 
orbit. 
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Figure 20. Generation of required control torques; z axis 
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Figure 21. Generation of required control torques; x axis 


The roll and yaw equations can be expressed in matrix form as 


Ta —KyBxBzAh |_| AC || 06) 78 
Ta, —KyByBzAh | | BD || y(s) ” 


where 


B= [(lez —lyy + Lax) @ on —h + K2BxBz]8 + [zz@ on + Ki BE] (79) 
C=—-[Ux - ly + Izz)@ on +h— K2B,Bz|s — [xx® on + K3B?] 


Using methods introduced in Reference 10, the nutation and orbital time constants of the 


system can be estimated by 


4hlyxlz2 


Tt ——— ee ee go 
7 BEA ext]ez)hK—laxle(4K3+K 1)] 
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4hwo 
i (80) 


4hoo 
UO) ee 
3mBo(K ,+2K3) 


where B, is M/R* of Equation (26) and the time constants are expressed in orbits. The 
values of K, and K; were selected to saturate the torque rods for an angular error of one 


degree. This leads to a selection of the control gains as 


K, =- 9.555 x 10’ Amp-m’/Tesla 
K, =- 4.777 x 10° Amp-m’ -sec/Tesla (81) 
K, = - 2.389 x 10’ Amp-m’/Tesla 


which results in time constant values of 
t, = 13.711 seconds 
tT, = 6.848 minutes (82) 


tT, = 5.326 minutes 


These time constants only apply under conditions of no torque rod saturation and 
continuous control (i.e. no eclipse shut off of attitude sensing data). Since the gains are 
sized to create torque rod saturation above one degree pointing error, the simulated time 
constants will be slightly longer. Of note, K, AND K, are destabilizing for nutation. 


This can be managed by ensuring that the following relationship is satisfied 
(4K;+K,) < [((41,.+1,)hKo\/1L, (83) 


Figures 22 and 23 show simulation results which can be used to estimate system 


time constants. Figure 22 shows nutation and Figure 23 shows the orbital time constants. 
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Using these figures to estimate the time constants yields a nutation time constant of 
approximately 21 seconds, and an orbital time constant of approximately 6.5 minutes. 


These estimations are close to the calculated values. 


pointing error in degrees vs time; sun-oriented contro! 
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Figure 22. Nutation time constant response 


Figures 24 through 34 show simulation results for sun-oriented control in different 
orbit configurations (1:30, 3:00, 4:30 and 6:00 orbits). Figures 24 through 27 do not 
include sensor noise inputs or cross product of inertia effects. Figures 28 through 31 do 
include these effects with a sensor noise gain of 0.07 degrees which is much higher than 
those typically found in commercial sensors. Figure 28 and Figures 32 through 34 
include results for a 1:30 orbit with noise and cross product of inertias at different Earth 
positions around the sun. Figure 28 simulates Summer Solstice, Figure 32 simulates 
Vernal Equinox, Figure 33 simulates Autumnal Equinox, and Figure 34 simulates Winter 
Solstice. All eleven simulations use the following initial conditions 


Oo 


6, = 2 degree Oyo = Dono radians / second (see Equation (53)) 


1 degree ®,) = 0.00075 radians / second 


Wy = —0.5 degree @,) = —0.00003 radians / second 
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Figure 24. Sun oriented control; 1:30 orbit; no noise input 


Figure 23. Orbital time constant response 
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Figure 25. Sun oriented control; 
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Figure 26. Sun oriented control; 
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Figure 27. Sun oriented control; 6:00 orbit; no noise input 


Figure 28. Sun oriented control; 
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Figure 29. Sun oriented control; 
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Figure 30. Sun oriented control; 
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pointing error in degrees vs time; sun-oriented control 
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Figure 31. Sun oriented control; 6:00 orbit; with noise input 
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Figure 32. Sun oriented control; with noise input. Vernal Equinox 
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pointing error in degrees vs time; sun-oriented control 
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Figure 33. Sun oriented control; with noise input. Autumnal Equinox 
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Figure 34. Sun oriented control; with noise input. Winter Solstice 
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V. CONCLUSIONS AND RECOMMENDATIONS 
A. CONCLUSIONS 


Equations of motion for a sun-oriented inertial reference frame were developed. 
Transformation to a body-oriented reference frame was accomplished using a 3-1-2 Euler 
axis rotation sequence with small angle approximations. 

Magnetic attitude control laws for roll and yaw control were developed by 
modifying control laws which have been used in the past for Earth-oriented operation. 
Pitch control was accomplished by varying momentum wheel speed using a proportional / 
derivative feedback controller. The Earth-orientation control laws work well for this 
sun-oriented control application. 

| Simulations were performed to test capture control of the satellite once it is 
separated from the launch vehicle. The satellite was separated wtih and without the 
momentum wheel being spun up to nominal operating speed to test system performance. 
The control law is capable of removing all rotational kinetic energy and aligning the 
momentum vector along the orbit normal. The possibility exists for the momentum vector 
to be aligned away from the sun (180 degree misdirection) when the satellite is separated 
before the momentum wheel 1s spun up. 

Reorientation from Earth-oriented / system check-out operations to a sun-oriented 
mode was performed using the same magnetic attitude control law as used for the 
sun-oriented roll and yaw control. The control law does not break down for low Beta 
angle orbit applications. 

Sun-oriented control maintains a pointing accuracy better than 0.2 degrees for the 
worst case, 1:30 orbit configuration. Given this, momentum wheel size could probably be 
cut in half to maintain the desired pointing accuracy within 0.5 degrees. This would 
reduce weight and cost while also reducing reorientation times from the Earth-oriented to 
the sun-oriented mode. A second Earth sensor is not required to maintain pointing 


accuracy during eclipse in the current configuration. 
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Noise and cross product of inertia inputs do not have a noticeable affect on 
pointing accuracy output. 
B. RECOMMENDATIONS 

Future study could be done which implements a Kalman Filter for noise 
attenuation in the sensors. The filter could also be used to estimate pitch rate from 
magnetometer data when Earth sensor data is unavailable. This would avoid the 
undamped oscillations experienced without pitch rate feedback. 

Study could be performed using a 6:00 orbit (continuous sun sensor input) without 
a momentum bias to test the ability of the torque rods to control the spacecraft by 
themselves. Torque rod size could also be modified to see how small the torque rods can 
be before they are unable to control the satellite. 

The current SIMULINK program could be modified in order to speed up 
simulation times by removing several of the function blocks and constructing “hard wired" 


representations. 
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APPENDIX: SIMULATION BLOCKS AND MATLAB CODE 


MATLAB/SIMULINK block diagrams are included here for the program used to 
simulate all satellite control sequences. Computer code for MATLAB control programs 
are also listed. The main SIMULINK program block is shown in Figure Al. The beta 
angle and Earth to Sun angle gernorater block creates beta angle and Earth to sun angles 
in accordance with Chapter II equations. The B-field block generates the Earth's 
magnetic dipole field in Earth-oriented orbital coordinates. This field is then transformed 
to body coordinates by performing the two coordinate transformations; 3-axis orbit to sun 
frame and 3-axis sun to body frame. The magnetometer block acts as the magnetic field 
sensor and time derivative estimator for the rate of change of the magnetic field with 
respect to the spacecraft body. These outputs are then fed to the Magnetic torque rod 
control along with change in wheel momentum (Ah), spacecraft attitude and spacecraft 
attitude rate. Magnetic moments are then created and fed to the torque rod block where 
saturation limiting is performed. The saturation limited magnetic moments are then used 
by the cross product block M x B to calculate magnetic control torque which 1s the cross 
product between the magnetic moments and the magnetic field expressed in body 
coordinates. These torques along with the disturbance torques are sent to a torque 
summer which feeds the overall external torque input to the spacecraft dynamics block 
where the body rotational rates (body rates) are calculated. The momentum from the 
momentum wheel is also input here to satisfy the equations of Chapter [V. The body 
rates are then input to the kinematics section to calculate the Euler angle rates which are 
integrated to form the Euler angles $, © and y. The attitude sensor block passes through 
attitude information with sensor noise input when sensor data is available. At all other 
times, is passes an error signal of zero to simulate the loss of attitude sensing information. 
Individual output blocks for each variable to be plotted are included to pass data to the 
MALAB workspace for subsequent processing. Each individual block is displayed on the 


following pages. 
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Figure A2. Earth's magnetic field in Earth-oriented coordinates 
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Figure A3. Beta and Earth to sun angle generation 
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Figure A4. Beta angle generation 
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Figure A5. Coordinate transformation from Earth orientation to Sun orientation 
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Figure A6. Coordinate transformation from Sun orientation to body orientation 
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Figure A7. Magnetometer sensor 
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Figure A9. Sun-oriented magnetic torque rod control 
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Figure Al5. Momentum wheel control 


a 
-_, 
Pitch Wheel V heel 


Vneel Wheel Command . 
Command DeMux Saturation Momentum 


| la 
al 


Gain2 Pitch Wheel Vheel 
torque Torque 
Saturation | 


wheeltorque 


To Workspace 





Figure Al6. Momentum wheel saturation controls 


73 














ce 
Sum2_ Integrator2 omega x °megas 
omega y 


omega z 





Momentum 


total momentum tolal momentum _ 


Figure Al7. Spacecraft Dynamics 


Integrator] 


Integrator2 


Integrator3 


omegaUnew 


o Workspace! ¢ 
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Figure A20. Gravity gradient torque generator 
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Figure A22. Generation of satellite velocity unit vector in body coordinates 
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Figure A23. Computation for Aero disturbance torque for each area 
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MATLAB CODE 


%this program is the controller for a simulink file projcon3 which will simulate 


“the equations of the attitude determination and control system 

%The program will set up several variables and initial conditions and will 

“%also control data display 

VL WLWLVAWAVASLYAVATAVAVANAVANATASANATATATASANANANANASANASANATASANA VATA TAVANATALATANAINA' 
printcon=0; 

control=1; “set to | to call stmulink program 
mode=0; 

lawcontrol=1; 

if control==1, 

clear 

control=1; 

lawcontrol=1; 

mode=0; 

printcon=0; 

end 

orient tall 

MAHMVMVMVMAMVMMMMMMM%%%%Y%%%V%0%0% 

%variable, constants and initial conditions 


maxrows=[80000,8]; %sets max # of data points 
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stepsizemin=.09; 
stepsizemax=.09; 

starttime=0; 

stoptime=10*5760; 

eclipsecon=0; 

rearth=6378.135; 

alt=550; 
earthrotvel=(360/24)*(1/3600)*(pi/180); 
wo=sart(398601.2/(rearthtalt)*3); 
rad/sec 

omega0=wo; 

inc=97.59306* pi/ 180; 
nu0=wo*starttime*p1/1 80; 

node 

epsilon=23.442* p1/180; 

% 

ieaneteis of inertia 

ixx=8.363; 

iyy=9.284; 

1zZ=4.012; 


ixy=-.111; 
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%set step size for simulation 


%1 means no eclipse 
Yradius of Earth in km 
%satellite altitude in km 
%spin rate of earth in rad/sec 


%orbital rotational velocity in 


%orbit inclination in degrees 


%angle around orbit from ascending 


%earth tilt wrt ecliptic 


%kg-m“2 





ixz=.013; 

1yz=.048; 

1i=[1XX -1Xy -1XZ}-Ixy Llyy -1YZ}-1XZ -1yz 1zz]; 

invI=inv(i1); “inverse of inertia matrix for dyn 

% 

“beta angle generation 

uoffset=88; %set inital Earth position around sun 

betaw0=(uoffset+22.5)*pi/180; 

betawdot=-9.9639*cos(inc)*(1/24)*(1/3600)*(pi/180)* (rearth/(rearth+alt))*3.5; 

betau0=(uoffset+0)*pi/1 80; 

betaudot=(360/365.242)*(1/24)*(1/3600)*(pi/180); 

beta0=asin(cos(inc)*sin(betau0)*sin(epsilon)+sin(inc)* sin(betaw0)*cos(betau0)... 
-cos(epsilon)*cos(betaw0)*sin(betau0)* sin(inc)); 

% 

%solve for nusun for earth to sun angle computations 

A=acos(-cos(pi-inc)*cos(epsilon)+sin(pt-inc)*sin(epsilon)*cos(betaw0)); 

phro=asin(sin(pi-inc)*sin(betaw0)/(sin(A)+.0000000001))-betau0; 

B=asin(sin(epsilon)*sin(betaw0)/(sin(A)+.00000000001)); 

bbb=tan(beta0); 

ccc=tan((pi/2)-A); 


eee=asin(bbb*ccc); 
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nusun=B-(eee); 

nuesQ=nu0+nusun; 

% 

%solve for theta earth to sun at time 0 
thetaesO=atan2(tan(nues0),sin(beta0)); 
AAAS ES ALSO SOLANA TA SATA NASA NANANASATANATAVASANATASASASANARASANASATASASASANA' 
“disturbance torques 

%center of mass listings in meters 
cex=.01*.0254; 

cgy=-.23*.0254; 

cgz=-17.84*.0254; 

% 

%disturbance torque control 

“aero cutoffs 

kal=1; 

ka2=1; 

ka3=1; 

ka4=1; 

ka5=1; 

ka6=1; 


ka7=1; 
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Yaxis cutoff by type of disturbance torque 

gexcon=1; 

solarxcon=1; 

aeroxcon=1; 

ggycon=1; 

solarycon=1; 

aeroycon=1; 

gezcon=1; 

solarzcon=1; 

aerozcon=1; 

% 

“axis torque cutoff 

distxcon=1; 

distycon=1; 

distzcon=1; 
UHVMVMM%MMMMMMMMMMMMMMMMHMNMMMMMNVMV%%/V%VVV%0%0%0/0 
“sensor noise 

earthsenscon=0*(.07*pi/ 180); %set nonzero for noise input 
rollsunsenscon=0*(.07*pi/180); 

yawsunsenscon=0*(.07* pi/1 80); 


% 


82 











%b-field generation initial conditions 

alfaQ=nu0; 

mu0=0; 

% 

%wheel control gains and momentum bias 

taup=10; 

Kp=3.6e-1; 

wheelbias=9.75; %wheel inertia is .0310352 kg-m“2 
wheelh0=0; 

% 

%torque rod control 

K1=(-5*1.911e7); 

K2=-5*9.554e7; 

K3=(K1)/4; 

Ku=-1*(1.538e6); 

AAAS VANA NASAL NA TATA SA NASATATATATANASATASASASASASASASASASASATATASASASA 
%kinematics initial conditions 

phi0=(1*pi/180); 

theta0=(2* pi/180); 

ksi0=(-.5* pi/180); 


WWW VPALWYVAVWDVV/VAWLVALLVAVVA%% 
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“spacecraft dynamics initial conditions 

wx0=0; 

wy0=((cos(thetaes0)“2)*(-wo)/((cos(nues0)*2)*sin(beta0))); 

wz0=0; 

% 

Hx0=1xx*wx0; 

Hy0=1yy*wy0-9.75; 

Hz0=1zz* wz0; 

VINA NANA NA NASA NASANAVANASATANASARASATANASANARAVANASATANATARASATASAVANASATASANASAVARLYAS! 
if control==1, 

% 

“call simulink program 

finalpr3 

pause 

end 

% 

wheelspeed=wheelmomentum*60/(2* pi*.03 10352); 
“percentspeed=100* abs(abs(wheelspeed)-3000)/300; 
Ypercentwheeltorque=100* ab s(hdotwheel)/ 025; 
Y%percentdeltamomentum=100*abs(wheeldeltamomentum/.9749996); 


Yo 


84 





Y%phisqrd=phi.* phi; 

%thetasqrd=theta.* theta; 

%ksisqrd=ksi.*ks1; 
pointerror=sqrt((phi.*phi)+(ksi.*ksi)); 
%%H%H%™%%H%H%%%%%M%M%M%%%%%M%%%%%%%%%%%%%%%VV%V%V0%0'/0/0 
Y%plot results 

clg 

if control==101, 
subplot(311),plot(time/60,bxorb) 

title(bfild in Earth-oriented frame’) 
xlabel(‘time in minutes’) 

ylabel(‘bfield strength x in Earth-orbit frame’) 
grid 

subplot(312),plot(time/60,byorb) 

xlabel(‘time in minutes’) 

ylabel(‘bfield strength y in Earth-orbit frame’) 
grid 

subplot(313),plot(time/60,bzorb) 

xlabel(‘time in minutes’) 

ylabel('bfield strength z in Earth-orbit frame’) 


grid 
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pause 

if printcon==1, 

print 

end 

% 

subplot(311),plot(time/60,bxsun) 

title(‘bfield in Sun-oriented frame’) 
xlabel(‘time in minutes’) 

ylabel(‘bfield strength x in sun-oriented frame’) 
grid 

subplot(312),plot(time/60,bysun) 

xlabel(‘time in minutes’) 

ylabel(‘bfield strength y in sun-oriented frame’) 
grid 

subplot(3 13), plot(time/60,bzsun) 

xlabel(‘time in minutes’) 

ylabel(‘bfield strength z in sun-oriented frame’) 
grid 

pause 

if printcon==1, 


print 
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end 

% 

subplot(311),plot(time/60,bxbody) 
title(bfield 1n body frame’) 

xlabel(‘time in minutes’) 

ylabel(‘bfield strength x in body frame’) 
grid 

subplot(3 12),plot(time/60,bybody) 
xlabel(‘time in minutes’) 

ylabel(‘bfield strength y in body frame’) 
grid 

subplot(313),plot(time/60,bzbody) 
xlabel(‘time in minutes’) 

ylabel(‘bfield strength z in body frame’) 
grid 

pause 

if printcon==1, 

print 

end 


end 


87 








LLU /AVW/VV/VALVU/ V/V/V ANNA Y0%% 
subplot(311),plot(time/60,02sphi* 180/pi) 
title(‘orb to sun angles vs time’) 
xlabel(‘time in minutes’) 

ylabel(‘o2s phi’) 

grid 

subplot(312),plot(time/60,02stheta* 1 80/p1) 
xlabel(‘time in minutes’) 

ylabel('o2s theta’) 

grid 

subplot(3 13),plot(time/60,02sksi* 180/p1) 
xlabel(‘time in minutes’) 

ylabel('o2s ksi') 

erid 

pause 

if printcon==1, 

print 

end 

% 

clg 


subplot(311),plot(time/60,02sthetadot) 


88 





xlabel(‘time in minutes’) 

ylabel('o2s thetadot’) 

grid 

subplot(312),plot(time/60,betaangle* 1 80/pi) 
xlabel(‘time in minutes’) 

ylabel('betaangle’) 

grid 

subplot(3 13),plot(time/60,woddotfunct) 
xlabel(‘time in minutes’) 
ylabel(‘omegaOnew dot rad/sec’2') 

grid 

pause 

if printcon—1, 

print 

end 

clg 

% 

if control==101 
subplot(311),plot(time/60,velunitvect(:,1)) 
title(‘velocity unit vector components in body frame’) 


xlabel(‘time in minutes’) 
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ylabel(‘velocity in x body') 

grid 

subplot(3 12),plot(time/60,velunitvect(:,2)) 
xlabel('time in minutes’) 

ylabel(‘velocity in y body’) 

grid 
subplot(313),plot(time/60,velunitvect(:,3)) 
xlabel(‘time in minutes’) 

ylabel(‘velocity in z body’) 

grid 

pause 

if printcon==1, 

print 

end 

end 

WL WANA WANA NASANLSASATASANASASANANATASATALAVASATASAVANATASLTASLTANAVASATASANASANASANALA 
subplot(311),plot(time/60,phi* 1 80/p1) 
title(euler angle in degrees’) 

xlabel(‘time 1n minutes’) 

ylabel(‘roll angle in degrees’) 


grid 
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subplot(312),plot(time/60,theta* 180/p1) 
xlabel(‘time in minutes’) 

ylabel(‘pitch angle in degrees’) 

grid 

subplot(313),plot(time/60,ksi* 180/pi) 
xlabel(‘time in minutes’) 

ylabel('yaw angle in degrees’) 

grid 

pause 

if printcon==1, 

print 

end 

% 

subplot(311),plot(time/60,phidot* 180/p1) 
title(‘euler angle rates in degrees per second’) 
xlabel(‘time in minutes’) 

ylabel(‘roll angle rate in deg/sec’) 

grid 

subplot(312),plot(time/60,thetadot* 180/pi) 
xlabel(‘time in minutes’) 


ylabel(‘pitch angle rate in degrees’) 
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grid 

subplot(3 13), plot(time/60,ksidot* 180/p1) 
xlabel(‘time in minutes’) 
ylabel('yaw angle rate in degrees’) 
grid 

pause 

if printcon==1, 

print 

end 

% 

if control==101, 
subplot(311),plot(time/60,ggx) 
title(‘grav grad torques vs time’) 
xlabel(‘time 1n minutes’) 
ylabel('grav grad x’) 

grid 
subplot(312),plot(time/60,ggy) 
xlabel(‘time in minutes') 
ylabel('grav grad y’) 

grid 


subplot(313),plot(time/60,ggz) 


92 








xlabel(‘time in minutes’) 
ylabel(‘grav grad z’) 

grid 

pause 

if printcon==1, 

print 

end 

%o 
subplot(311),plot(time/60,aerox) 
title(‘aero torques vs time’) 
xlabel(‘time in minutes’) 
ylabel(‘aero x’) 

orid 
subplot(312),plot(time/60,aeroy) 
xlabel(‘time in minutes') 
ylabel(‘aero y’') 

grid 

subplot(3 13),plot(time/60,aeroz) 
xlabel(‘time in minutes’) 
ylabel(‘aero z’) 


erid 
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pause 
if printcon==1, 

print 

end 

end 

% 

clg 

subplot(411),plot(time/60,wheelmomentum) 
title(‘momentum wheel characteristics vs time’) 
xlabel(‘time in minutes’) 

ylabel(‘wheel momentum’) 

grid 

subplot(412),plot(time/60,hdotwheel) 
xlabel(‘time in minutes’) 

ylabel(‘h dot[ wheel torque]’) 

erid 
subplot(413),plot(time/60,wheeldeltamomentum) 
xlabel(‘time in minutes’) 

ylabel(‘delta momentum’) 

grid 


subplot(414),plot(time/60,wheelspeed) 
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xlabel(‘time in minutes’) 

ylabel('wheel speed in rpm’) 

grid 

pause 

if printcon==1, 

print 

end 

% 

clg 

subplot(3 11),plot(time/60,magmomentsout(:,1)) 
title(‘magnetic torque rod moments vs time’) 
xlabel(‘time 1n minutes’) 

ylabel(‘actual mag moment x’) 

grid 
subplot(312),plot(time/60,magmomentsout(:,2)) 
xlabel(‘time in minutes’) 

ylabel(‘actual mag moment y’) 

grid 

subplot(3 13),plot(time/60,magmomentsout(:,3)) 
xlabel(‘time 1n minutes’) 


ylabel(‘actual mag moment z’) 
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grid 

pause 

if printcon==1, 

print 

end 

% 
subplot(311),plot(time/60,magtorquex) 
title(magnetic torque rod torques vs time’) 
xlabel(‘time in minutes’) 

ylabel('mag torque x’) 

grid 

subplot(3 12),plot(time/60,magtorquey) 
xlabel(‘time in minutes') 

ylabel(‘mag torque y’) 

grid 
subplot(313),plot(time/60,magtorquez) 
xlabel(‘time in minutes’) 

ylabel('mag torque z') 

grid 

pause 


if printcon==1, 
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print 

end 

% 
subplot(311),plot(time/60,totaltorquein(:,1)) 
title(‘total torques each axis vs time’) 
xlabel(‘time in minutes’) 

ylabel(‘total torque x in N-m’) 

grid 

subplot(3 12),plot(time/60,totaltorquein(:,2)) 
xlabel('time in minutes’) 

ylabel('total torque y in N-m’) 

grid 

subplot(3 13),plot(time/60,totaltorquein(:,3)) 
xlabel(‘time in minutes’) 

ylabel(‘total torque z in N-m’) 

grid 

pause 

if printcon==1, 

print 

end 


% 
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subplot(311),plot(time/60,disttorquex) 
title(‘dist torques each axis Vs time’) 
xlabel(‘time in minutes’) 

ylabel(‘dist torque x in N-m’) 

grid | 

subplot(3 12),plot(time/60,disttorquey) 
xlabel(‘time in minutes’) 

ylabel('dist torque y in N-m') 

grid | 
subplot(3 13),plot(time/60,disttorquez) 
xlabel(‘time in minutes’) 

ylabel(‘dist torque z in N-m’') 

grid 

pause 

if printcon==1, 

print 

end 

% 
subplot(311),plot(time/60,bfldsensedbody(:,1)) 
title(‘sensed b-field vs time’) 


xlabel(‘time in minutes’) 
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ylabel(‘sensed b-field x dir’) 

grid 
subplot(312),plot(time/60,bfldsensedbody(:,2)) 
xlabel(‘time in minutes’) 

ylabel('sensed b-field y dir’) 

grid 
subplot(313),plot(time/60,bfldsensedbody(:,3)) 
xlabel(‘time in minutes’) 

ylabel('sensed b-field z dir’) 

grid 

pause 

if printcon==1, 

= 

end 

% 

%end 

% 

clg 

subplot(3 11),plot(time/60,omegaout(:,1)) 
title(‘body rates’) 


xlabel(‘time in minutes’) 
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ylabel(‘omega x’) 

grid 
subplot(312),plot(time/60,omegaout(:,2)) 
xlabel('time in minutes’) 

ylabel(‘omega y') 

grid 

subplot(3 13), plot(time/60,omegaout(:,3)) 
xlabel(‘time in minutes’) 

ylabel(‘omega z’) 

grid 

pause 

| if printcon==1, 

print 

end 

% 

clg 
subplot(311),plot(time/60,totalmomentum(:,1)) 
title(‘total momentum’) 

xlabel(‘time in minutes’) 

ylabel(‘total momentum x’) 


grid 
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subplot(3 12),plot(time/60,totalmomentum(:,2)) 
xlabel(‘time in minutes’) 

ylabel(‘total momentum y’) 

grid 

subplot(3 13),plot(time/60,totalmomentum(:,3)) 
xlabel(‘time in minutes’) 

ylabel(‘total momentum z’) 

grid 

pause 

if printcon==1, 

print 

end 

% 

cle 

plot(time/60,pointerror* 1 80/pi,'b') 
title(‘pointing error in degrees vs time; sun-oriented control’) 
xlabel(‘time 1n minutes’) 

ylabel(‘pointing error to sun’) 

grid 

hold 


plot(time/60,switchout,'b--') 
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pause 
if printcon==1, 
print 


end 
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